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Abstract 

This paper addresses the fine-scale axisymmetric structure exhibited in Saturn's 
A and B-rings. We aim to explain both the periodic microstructure on 150-220m, 
revealed by the Cassini UVIS and RSS instruments, and the irregular variations 
in brightness on 1-lOkm, reported by the Cassini ISS. We propose that the former 
structures correspond to the peaks and troughs of the nonlinear wavetrains that 
form naturally in a viscously overstable disk. The latter variations on longer scales 
may correspond to modulations and defects in the wavetrains' amplitudes and wave- 
length. We explore these ideas using a simple hydrodynamical model which captures 
the correct qualitative behaviour of a disk of inelastically colliding particles, while 
also permitting us to make progress with analytic and semi-analytic techniques. 
Specifically, we calculate a family of travelling nonlinear density waves and de- 
termine their stability properties. Detailed numerical simulations that confirm our 
basic results will appear in a following paper. 
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1 Introduction 



Saturn's A and B-rings sport an abundance of irregular radial structure which, 
though aesthetically pleasing, presents something of a puzzle to the theoreti- 
cian. The instruments aboard the Cassini space probe show that these pat- 
terns manifest on a vast range of length-scales and can take quite different 
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forms. For instance, there exist quasi-periodic microstructure on scales of 0.1 
km (Colwcll ct al. 2007; Thomson ct al. 2007), discontinuous and irregular 
striations on the 1-10 km intermediate scale, and much broader 100 km undu- 
lations (Porco at al. 2005). In addition to the difficulties involved in tackhng 
these three orders of magnitude, there are the formidable modelling questions 
posed by a cold disk of densely-packed, inelastic, and infrequently colliding 
particles (Stewart ct al. 1984, Araki and Tremaine 1986, Salo 1991, Hameen- 
Antilla and Salo 1993, Schmidt et al. 2001, Latter and Ogilvie 2008). This 
paper will focus on only a subset of these phenomena, the smaller-scale vari- 
ations, and will not linger especially on the modelling issues. Specifically, it 
investigates how the quasi-periodic microstructure relates to the viscous over- 
stability, on one hand, and to structure formation on the intermediate scales, 
on the other. The variations on the much longer 100km scale are not examined, 
and we suspect that they have their origin in a different mechanism entirely, 
perhaps baUistic transport (Durisen 1995). 

Our starting point is the viscous overstability, which is now regarded as a key 
player in the short scale radial dynamics of Saturn's rings (Schmit and Tschar- 
nuter 1995, Schmidt et al. 2001, Spahn and Schmidt 2006, Latter and Ogilvie 
2008). The viscous overstability is an axisymmetric oscillatory instability that 
affiicts the homogeneous state of Keplerian shear. Growing modes rely on the 
alliance of the fluid disk's inertial-acoustic oscillations with the disk's stress 
oscillations: variations in the stress extract energy from the Keplerian shear 
and inject it into the inertial-acoustic wave; but the increased motion this 
induces magnifies the stress oscillation itself which can extract even more en- 
ergy, and so the process runs away. The feedback loop requires (a) the stress 
to efficiently remove energy from the Keplerian flow and (b) the two oscilla- 
tions to communicate effectively, in particular for them to be in phase. The 
first condition is tied to the stress's sensitivity to surface density. The second 
condition is often violated in dilute rings and turbulent disks, where the stress 
can lag behind the epicycles and the overstability fails to work (Ogilvie 2001, 
Latter and Ogilvie 2006a). 

The instability's linear theory has been well established by a variety of the- 
oretical approaches: hydrodynamics (Schmit and Tscharnuter 1995, Schmidt 
et al. 2001), A^-body simulations (Salo 2001, Salo et al. 2001), and kinetics 
(Latter and Ogilvie 2006a, 2008). However, its nonlinear theory has received 
surprisingly little attention. Two hydrodynamical studies exist, a large-scale 
nonlinear simulation of a ring annulus in ID (Schmit and Tscharnuter 1999) 
and a weakly nonhnear analysis (Schmidt and Salo 2003). The simulations 
show that the nonlinear evolution of an overstable disk is characterised by sig- 
nificant disorder. In contrast, the latter study suggests that an overstable ring 
may exhibit simple coherent structures which take the form of travelling waves. 
Because Schmit and Tscharnuter impose refiecting boundary conditions, and 
hence break translational symmetry, such coherent structures may have been 
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difficult to observe in their simulations. Certainly, more simulations need to be 
undertaken and the influence of the boundary conditions better understood. 
Concurrently, a fully nonlinear theory extending the work of Schmidt and 
Sale is required to fully estabhsh the existence of the nonhnear solutions. The 
former project we present elsewhere, the latter we present here. 

First, we demonstrate that steady nonlinear travelling wavetrains are exact so- 
lutions to the governing nonlinear equations of a viscously overstable disk, and 
second, that these solutions are the loci of a rich secondary set of dynamics 
which generate irregular variations in the waves' amplitude and wavenum- 
ber. We propose that viscously overstable regions in the A and B-rings sup- 
port a bed of travelling wavetrain solutions and these structures correspond 
to the quasi-periodic variations registered by Cassini's UVIS and RSS instru- 
ments (the radial 'microstructure'). Modulations and defects in the amplitudes 
and wavelengths of these wavetrains may yield different optical properties 
which can be traced by the Cassini cameras. Consequently, we hypothesise 
that these larger-scale variations are associated with the intermediate 1-10 
km structures observed (Porco et al. 2005). Our hypotheses are investigated 
with a one-dimensional, isothermal fluid model endowed with a Navier-Stokes 
stress. Though simplistic, it should predict qualitatively correct behaviour, 
while permitting the problem to be attacked analytically or semi-analytically. 
Importantly, self-gravity is omitted throughout, but we investigate its role in 
a following paper. This is mainly for readability as the nonself-gravitating 
dynamics is sufficiently complex on its own. 

A study of nonlinear waves connects naturally to the extensive field of wave 
propagation in thin astrophysical disks (Kato 2000), and in particular to the 
launching of spiral density waves (Goldreich and Tremaine 1978a, 1979, 1980, 
Borderies et al. 1983, 1986, Shu et al. 1985a, 1985b). In the latter studies the 
emphasis has been on the launching of spiral waves by an external potential, as 
might issue from a moon, and the damping of such waves by a viscous stress. 
But the viscous stress can be an engine of growth also, as in the development of 
global eccentric modes and narrow ringlets (Borderies et al. 1985, Papaloizou 
and Lin 1988, Lyubarskij et al. 1994, Longaretti and Rappaport 1995, Ogilvie 
2001) and the axisymmetric viscous overstabihty itself (Kato 1978, Schmit and 
Tscharnuter 1995). In this paper it will be shown that the viscous forces in a 
planetary ring can balance energy dissipation and injection and thus sustain 
steady wave structures. 

A study of nonlinear wavetrains also connects to the perhaps more massive 
field of nonlinear waves in general media and ID reaction-diffusion systems 
especially, the model equation of which is the complex Ginzburg-Landau equa- 
tion (Aranson and Kramer 2002). The dynamics of the complex Ginzburg- 
Landau equation is remarkably rich and may provide a template for the study 
of viscous overstabihty in Saturn's rings. As in a fluid disk, the equation ad- 
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mits a trivial homogeneous sohition that is susceptible to an oscillatory linear 
instability via a Hopf bifurcation (the analogue of the viscous overstability); 
it also supports both stable and unstable steady nonlinear wavetrains. Varia- 
tions upon the wavetrains exhibit a wide variety of behaviours ranging from 
smooth modulations to abrupt jumps in wavenumber and amplitude (sources 
and shocks), as well as small-scale chaotic fluctuations (Bcrnoff 1988, Shraiman 
et al. 1992, Chate 1994, Aranson and Kramer 2002). In particular, the sources 
and shocks can partition the radial domain into a one-dimensional cellular 
structure, in which each 'cell' is characterised by a different wavetrain (Chate 
1993, Popp et al. 1994, Ipsen and Hecke 2001). The latter behaviour seems 
especially pertinent to ring structure on intermediate scales, which Cassini's 
cameras show as a pattern of irregularly spaced, disjunct bands (Porco et 
al. 2005). 

A summary of the paper is as follows. First, space will be lent to a discussion 
of the validity of the simple hydrodynamic model that we employ, the parame- 
ters it introduces, and the values these should take. In Section 3, which comes 
next, a brief summary of the linear stability of the homogeneous state of Kep- 
lerian shear will be presented. Though these results have appeared a number 
of times elsewhere (Schmit and Tscharnuter 1995, Schmidt et al. 2001), we 
include them for completeness. Section 4 will demonstrate the existence and 
investigate the properties of axisymmetric, steady, nonlinear, travelling wave- 
trains in a viscously overstable disk. The nonlinear wave profiles are calculated 
numerically via the solution of a nonlinear eigenvalue problem and analyti- 
cally in the limit of long wavelength. We find that an overstable disk supports 
a one-parameter family of wavetrain solutions, the members of which can be 
parametrised by their wavcmimbcr A. In Section 5 we establish the linear 
stabihty of the wavetrains and find that above a critical wavelength Agt the 
wavetrains are linearly stable. For model parameter values, associated with 
optical depths of 1-1.5, the critical wavelength is approximately 50H, where 
H is the scale height of the disk. The critical wavelength appears consistent 
with Cassini observations of periodic microstructure. Lastly, in Section 6, we 
discuss more fully the general nonlinear dynamics of an overstable disk. A 
rough argument is sketched explaining the upward cascade to longer scales 
observed in the simulations of Schmit and Tscharnuter (1999) and we sug- 
gest that, irrespective of self-gravity, this 'inverse cascade' should halt (or at 
least change qualitatively) once most power reaches a Icngthscale near Agt- A 
discussion follows which addresses what this saturated state may look like. 
Modulations and defects in the wavetrains' wavenumber and amplitude are 
touched on, and a multiple-scales analysis presented in which we show that 
sources and shocks in the wavetrains' phase may develop. Finally wc discuss 
the relevance of the complex Ginzburg-Landau equation. Section 7 presents 
our conclusions, in which we summarise the paper, and point out the issues 
which remain unresolved and necessary future work. 
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2 The model 



In order to bring out the important features of the ring's nonhnear evolution 
we employ a simple hydrodynamical model that captures the correct qualita- 
tive behaviour while not burdening the analysis or obscuring its interpreta- 
tion with mathematical complexity. The ring is assumed to be a vertically- 
averaged, non-self-gravitating, isothermal, ideal gas possessing a Newtonian 
viscous stress. In addition, the shearing and bulk viscosities depend on sur- 
face density as power laws, in accordance with Schmit and Tscharnuter (1995, 
1999) and Schmidt et al. (2001). We also work with the shearing sheet ap- 
proximation, which is a Cartesian representation of a 'patch' of disk orbiting 
the central planet with frequency fl, and in which x and y denote the radial 
and azimuthal dimensions respectively (see Goldreich and Lynden-Bell 1965). 

The governing equations are 

dta + d,{aui) = 0, (1) 
a{dtUi + UjdjUi + 2QeizjUj) = —adi^r — diP + dj^ij-, (2) 

where a, P, and Uij are the vertically integrated density, pressure, and vis- 
cous stress, and Ui is the planar fluid velocity. The tidal potential is $t — 
-3Q2xV2. 

The pressure comes from the ideal gas equation of state 

P = v',a, (3) 
where Vs is the isothermal sound speed. The viscous stress is given by 

Uij = iya{diUj + djUi) + {ub - lu)a{dkUk)Sij. (4) 
The kinematic shear and bulk viscosities are parametrised as 

where a, ai,, and (3 are dimensionless parameters, and cr* is the surface density 
of the homogeneous state of Keplerian shear. 



Throughout the paper we employ dimensions so that 

Q — 1, Vg — i, (J* — 1. 

One unit of time is therefore (27r)^^ times an orbital period and the unit of 
length is H = Vs/^, the scale height of the disk. The full set of governing 
parameters is now a, aj,, and (3. 
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2.1 Assumptions 



Of course, much can be said about each of the model assumptions. We shall say 
but a little. First, the adoption of vertical averaging limits us to radial scales 
much longer than the disk scale height, H. The shearing sheet on the other 
hand limits us to scales much shorter than the disk radius. Neither constraint 
is much of a problem because the lengthscales of the phenomena we examine 
sit well within this enormous range. More of an issue is the neglect of the 
vertical motions. For instance, the viscous overstability is usually accompanied 
by vertical 'breathing' or 'splashing' and may be significant when combined 
with nonisothermal behaviour. We, however, leave this issue open for the time 
being. 

A real planetary ring is not generally isothermal. Due to the relative infre- 
quency of collisions, the thermal time-scale is of the same order as the dynam- 
ical time-scale, as kinetic theories and iV-body simulations show (Goldreich 
and Tremaine 1978b, Stewart et al. 1984, Brahic 1977). That said, when col- 
lisions are a little more frequent — as in a dense ring which may support 
some tens of coUisions per orbit — the isothermal model affords an acceptable 
approximation (Salo et al. 2001, Latter and Ogilvie 2008). 

A dense ring is not an ideal gas, and is probably better suited to a polytropic 
equation of state which can better capture the effects of close packing and 
nonlocal pressure (Schmidt et al. 2001, Schmidt and Salo 2003, Latter and 
Ogilvie 2008). We persist, however, with the simpler model in the belief that 
the errors introduced do not alter the qualitative behaviour. 

Finally, it should be acknowledged that the viscous stress of a planetary ring 
behaves in a way not always captured by a simple Navier-Stokes stress. For 
instance, in a dilute ring, effects associated with nonlocality in time ('memory 
effects') are important because the translational stress relaxes on a time-scale 
comparable to the dynamical time-scale (Latter and Ogilvie 2006a). In con- 
trast, a dense ring possesses a stress dominated by the coUisional component, 
and, though local in time, will be a nonlinear function of the rate of strain 
(Latter and Ogilvie 2008). The Navier-Stokes model, however, offers a reason- 
able approximation of the dense system, particularly as we may 'tweak' the 
parameter (3 in order to roughly 'compensate' for the linearisation in strain. 

Self-gravity is excluded in this paper and will be examined separately in a 
subsequent article. It undoubtedly plays a role in both the axisymmetric and 

nonaxisymmetric dynamics of an overstable disk but the effects are compli- 
cated and deserve a special treatment. In addition, A^-body simulations show 
that nonaxisymmetric self-gravity wakes hinder the development of the linear 
modes (Schmidt et al. 2001), and it is most probable that they impact on the 
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nonlinear dynamics as well. This is an important issue but one that can only 
be resolved satisfactorily by three-dimensional simulations. 



2.2 Parameters 



Aside from these theoretical issues, we are encouraged to use the isothermal 
hydrodynamic model because it predicts qualitatively correct behaviour when 
compared with A^-body simulations. The linear growth rates of the overstable 
modes arc adequately approximated (Schmidt et al. 2001) as is the weakly 
nonlinear cvohition of their amplitudes (Schmidt and Salo 2003). For these 
reasons the hydrodynamic parameter values we use will be set equal to those 
computed in these studies, in particular, from Salo et al.'s (2001) A^-body 
simulations, which treated particles of radius 100 cm, undergoing collisions 
according to the Bridges et al. (1984) piecewise power law at a radius of 
100, 000 km. Full self-gravity was not included in the simulations, but its 
compression of the disk thickness was mimicked by increasing the vertical 
epicyclic frequency of the particles (an idea pioneered by Wisdom and 
Tremaine 1988). In Table 1 we reproduce some of the data of these runs 
for different optical depth r and with a vertical frequency enhancement of 

n^/n = 3.6. 

These values will serve us only as a guide. The viscous parameters are closely 
tied to the kinetic parameters of a particular simulation (such as particle size, 
elasticity law) but their complete functional dependence is far from under- 
stood. On the other hand, the kinetic parameters of a real particulate ring are 
not yet fully constrained. Thus we let a, ah and especially j3 take a variety 
of values. The j3 parameter is the quantity most sensitive to the background 
optical depth, as we can see from Table 1, and also to self-gravity. With re- 
spect to the last point, Schmidt and Salo (2003) have pubhshed another set 
of parameter values from a simulation in which il^/Jl — 2, not 3.6, and they 
find that j3 takes significantly smaller values in this case. On the other hand, 
the simplifying assumption of a Newtonian stress may require us to vary (3; 
and indeed, Schmidt and Salo (2003) infiate /3 by up to 13% in order to ob- 
tain agreement in their weakly nonlinear analysis. Finally, in a real disk, the 
existence of gravitational wakes (screened out by the A"-body simulations we 
mention) will militate against the development of the overstable modes, mean- 
ing that j3 may take a smaller 'effective' value. The actual situation, of course, 
is probably more complicated, not least by the additional effects of nonlocal 
'gravitational viscosity' (Daisaka et al. 2001). 
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T 


a 


Q!5 


R 


Vg/^ (m) 


yj.o 




1 ns 


u.u / 


9 47 


1.0 


0.357 


0.764 


1.15 


3.29 


1.5 


0.342 


0.681 


1.19 


4.42 


2.0 


0.322 


0.683 


1.55 


5.45 



Table 1 

Nondimensionalised hydrodynamical parameters at different optical depths r de- 
rived from AT-body simulations for a disk of 100 cm radius particles, at a radius of 
100, 000 km, undergoing collisions according to the Bridges et al. (1984) elasticity 
law, with vertical frequency enhancement of il^/fi = 3.6 (Schmidt et al. 2001, Salo 
et al. 2001). 



3 Linear stability of the homogeneous steady state 



This section presents a summary of the stability analysis of the homogeneous 
steady state of Keplerian shear. Though it has been thoroughly examined with 
various continuum models ( Lin and Bodenheimer 1981, Ward 1981, Stewart et 
al. 1984, Schmit and Tscharnutcr 1995, Spahn et al. 2000, Schmidt et al. 2001, 
Latter and Ogilvie 2006a, 2008) we repeat the analysis for completeness and 
to connect it to our nonlinear work in Sections 3 and 4. 

The equations (1) and (2) admit the following equilibrium: a = 1 and Uy ~ 
— (3/2) X, Ux = 0. We introduce a small perturbation so that a — 1 -\- a' and 
u = — (3/2)xej,-|-u'. Their hnearised equations read 

dt(j' + 9X = 0, (6) 

- 2u'y = -dj', (7) 
dtu'y + i< = -dxh', (8) 

with 

3 

f = a' - (ttfe + ^a)dxu'^, h' = -a{l + f3)a' - a d^Uy. (9) 

(For a full list of symbols see Table 2.) The first term in the expression for 
h! is the key to the disk's stabihty. It couples the variations in the viscous 
stress with the background shear and thus allows energy to be drawn from the 
shear and directed into the perturbation. We now represent the perturbation 
as a Fourier mode oc e*^^"*"** where k is the (real) wavenumber and s is the 
(complex) growth rate. The ansatz is substituted into Eqs (6)-(9) from which 
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Ibymbol 


Definition 


r 


Normal geometric optical depth 


a 


buriace density 


u, u' 


Total and perturbation velocities 




lidai potential 


/ 


Pressure' 


TT 
11 


Viscous stress 


J, g 


Dimensionless xx and a;y component of the pressure tensor 




Shearing and bulk kinematic viscosities 


a, ab 


Shearing and bulk 'alpha' parameters 


n 

p 


Exponent of density dependence of the kinematic viscosities 


Vs 


Sound speed 


\l 


Orbital frequency 


TT „. /O 

H = Vg/il 


Disk scale height 


k, s 


Wavenumber and growth rate of a linear disturbance 




Wavenumber and frequency of a nonlinear wavetrain 


Cp = Uj/fJ, 


Phase speed of a nonlinear wavetrain 


Cg = duj/dfi 


Group speed of a nonlinear wavetrain 


9 = fix — ut 


Phase variable of a nonlinear wavetrain 


A = Z7r//x 


Wavelength of a nonlinear wavetrain 


Ast 


Critical wavelength above which nonlinear wavetrains are linearly stable 


Z = [a,u^,Uy) 


Solution vector 




Dimensionless velocity perturbation variables 




Amplitude and phase variables for asymptotic long wavelength wavetrains 


r 



Small ordering parameter 


X = (){.}■ - Cyl). T = 5-1 


Long si")Mlial and l('ni|")oral \'arial)los 


e{x,T),K = ex 


Slowly varying phase and wavenumber of a nonlinear wavetrain 


Csh 


Speed of a travelling sliockfront 



Table 2 

Table of parameters, variables, and symbols 

one may derive the following dispersion relation, 



1 + /c^ + k'^a^^a + «{,) 



4 

3' 



3{l + (3) + k' 



(10) 



In the long wavelength limit, < /c <^ 1, two instabilities may be distin- 
guished: the viscous instability and the viscous overstability. The former pos- 
sesses the growth rate 

and is unstable if /? < —1. The instability will be extinguished on short wave- 
lengths and the critical k upon which marginal stability resides can be com- 
puted from \k\ — ^—3(1 -|- Wavelengths shorter than this value are vis- 
cously stable. Table 1 states that a dense particulate ring possesses a positive 
P for all T, and so the viscous mode will probably decay in real dense planetary 
rings. 

The viscous overstabihty, on the other hand, possesses growth rates 



±i 



1 + Ik^] + \ [3a(l + /3) - aa + aA k^ + 0(|A;|3), (11 



and is hence oscillatory, corresponding to either standing waves or right or left- 
going travelling waves. In the long wavelength regime, overstability emerges if 
the following criterion is satisfied: 

/'>ir^-?i. (12) 



3 V a 3 

(Schmit and Tscharnuter 1995). More generally, we have /3 > /3c, where /3c 
summarises the thermal properties of the system. A naive application of the 
condition (12) to the data of Table 1 indicates that dense planetary rings with 
optical depths of at least 1 are viscously overstable (but compare with the 
detailed kinetic results of Latter and Ogilvie 2008). Like the viscous instabihty 
sufficiently short scales will be stabihsed by pressure. The wavenumber of the 
marginal modes can be computed from the quartic equation 

a{a.}) -|- \oL){pH) + \(x)k'^ + (ctft + |q;)A;^ + [a.^, — \a — 3q;/3) = 0. (13) 



4 Exact nonlinear solutions: uniform wavetrains 



Jf..l Introduction 



The linear theory is easy to establish; a more challenging task is to determine 

how the system evolves once it enters the nonlinear regime. At first, small 
perturbations will grow independently and exponentially in the form of over- 
stable modes, but when their amplitudes become sufficiently large they will 
interact. At this point the evolution of the system becomes difficult to track 
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and, as with many hydrodynamical problems, we may be obliged to employ 
numerical simulations to describe its temporal and spatial variations. 

Only one fully nonlinear hydrodynamical study has so far been undertaken, 
that of Schmit and Tscharnuter (1999). They numerically evolved an isother- 
mal, one-dimensional annulus of gas unstable to the viscous overstability. The 
temporal extent of their simulations reached some 10,000 orbits and the radial 
extent some 4,000// (corresponding to about 40 km). These showed that the 
evolution of an overstable ring is disordered, but nevertheless manifests a char- 
acteristic spatial scale and a longer 'beating' pattern (though this 'beating' is 
perhaps an artefact of the finite size of the box). In addition, an overstable 
system exhibits an upward cascade of power to larger scales. If self-gravity 
is included, this process halts once a certain lengthscale is reached (equal to 
about 24i/ if a = = 0.2628 and /3 = 1.26). If self-gravity is omitted they 
report that the upward transfer of power does not cease: by 10,000 orbits the 
system has injected most of its power to a scale of some lOOH. 

An important feature of these simulations is their enforcement of reflecting 
boundaries. Though the domain is large, these boundary conditions will ulti- 
mately impede the formation of structures associated with translational sym- 
metry, such as travelling wave trains (at least globally). This is an impor- 
tant point because the weakly nonhnear analysis of Schmidt and Salo (2003) 
predicts that an overstable disk may support steady nonlinear waves, as do 
A'"-body simulations (Salo and Schmidt 2007). 

We take the view that the saturation of the overstability will indeed be ir- 
regular (as in Schmit and Tscharnuter 1999) but that these spatio-temporal 
variations will manifest upon a bed of coherent nonlinear wavetrain solutions 
(similar to those revealed by Schmidt and Salo 2003). The latter solutions we 
may regard as 'fixed points' in the disk's state space, i.e. the simplest non- 
trivial invariant solutions of the dynamics. The system's state trajectories we 
expect to wander around these points and in so doing exhibit much of their 
coherent structure. Similar phenomena characterise a number of diverse phys- 
ical systems such as reaction-diffusion systems, flame fronts, and ocean waves 
(Whitham 1974, Infeld and Rowlands 1990, Doelman et al. 2009), as well as 
higher dimensional turbulent systems, such as pipe Poiseuille flow and plane 
Couette flow (Gibson et al. 2008 and references therein). 

In this paper we make a start on investigating this idea. Our first task is 
to formally establish the existence of the axisymmetric nonlinear travelling 
wavetrains, the fixed points, from the fully nonlinear equations (l)-(2). This 
is what we undertake in this section. Once this is done we can probe their 
linear stability (Section 5) and then speculate on their role in the nonlinear 
dynamics generally (Section 6). 
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4-2 Governing equations 



The axisymmetric nonlinear perturbation equations are 



dta + u'^dxO + ad^u'^ = 0, (14) 

" a 

dtu'y + u'^d^u'y + -u'^^ a^^/i, (16) 



with 



/ = a-(«fe + |a)ai+/^aX, (17) 
h^-aa'+^{-l + d,u'y). (18) 

Wc propose that there exist wavetrain solutions with frequency uj and wavenum- 
ber II, and which move with the phase speed Cp = cu/ii. The waveprofiles are 
stationary in a frame moving with this speed, and so we let the density and 
velocity fields depend on the single phase variable: 

e^lix-cvt. (19) 

In addition, so as to ensure a uniform wavetrain, it is assumed that the density 
and velocity are 27r periodic in 6. In summary, we have 

Z={a,u'^,u'y), Z{x,t)^Zo{0), Zo(0) = Zo(27r). (20) 

Henceforth we drop the subscript 0. Substituting this ansatz into the continu- 
ity equation (14) supplies a first integral: 

o-{u'^ — Cp) = est. (21) 

The constant can be computed from angular momentum conservation: if we 
integrate the y-momentum equation over one period in 9 we obtain 



1 r2. 

2 Jo """^ 



dO = 0. 



Integration of (21) yields est = — Cp, which recommends a new dependent 
variable u — Cp — u'^. The surface density can then be neatly expressed: 

(7 = ^. (22) 
u 

Also we relabel v = u'y for notational uniformity. 

In order to compute the waveprofiles, u and v, we must solve the momen- 
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turn equation for the perturbations. This can be expressed as four ordinary 
differential equations: 




df du 2 



(23) 



dh dv Cp 



(24) 



de de 2n 



du / — cr 



(25) 



de /x(a6 + |Q)ai+/3' 
dv 1/3 ^ \ 



de IX \2 aa^+P ) ' 



(26) 



This set is subject to the periodic boundary conditions presented in (20). 

The system (22) — (26) is a nonhnear eigenvalue problem with input parame- 
ters a, ai,, /3, and fi and with eigenvalue co. Its solution can be accomplished 
numerically. These solutions we present in Section 4.5, but first we estab- 
lish some integral relations and some asymptotic results which elucidate more 
clearly the various processes and balances at work. 



4-3 Integral relations 

Two integral identities aid in the understanding of the solutions. The first es- 
tablishes the energy balance of the steady waves and shows that the important 
terms are those responsible for viscous dissipation and viscous overstability. 
The second provides an expression for the angular momentum transported by 
the waves. 

We derive two energy equations for the perturbations. The first governs the 
combined radial kinetic energy and thermal energy: 



where Uij is the viscous stress. The second governs the azimuthal kinetic 
energy: 



We now assume a steady wavetrain solution with wavelength A and integrate 
over one wavelength. The time derivatives and flux terms on the left sides of 
Eqs (2 7)- (28) vanish and we are left with the integrated source terms on the 
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right sides: 



r\ f-A 

[2au'^u'y - {d^u'^)Ii^:c\ dx = 0, [\<^K'^'y + {dxu'y)Ii^y\ dx = 0. 

These are combined so as to remove the 'Reynolds stress' term, and this 
furnishes us with the key balance: 

a { Qvd.u'y - [{v, + {d,u',f + Av{dxu'yf] } = 0. (29) 

In order to satisfy this equation, the first term, associated with viscous over- 
stability, must integrate to a positive value and balance the second viscous dis- 
sipation term which is always negative. The equation is an energy constraint 
that a uniform amplitude wavetrain must satisfy: kinetic energy removed by 
viscous dissipation must be exactly balanced by kinetic energy injected from 
linear viscous overstability. As such, Eq. (29) provides a means by which we 
can calculate the rate of strain (and nonlinear amplitude) of a wavetrain. 

An expression for the total angular momentum flux over a wavelength is easy 
to derive. We find the flux is equal to 

- / {au'^u' - U^y) dx = — {h-Cpv) dO. (30) 

A Jo ZTT Jo 

where we have used (23). In order to isolate the contribution from the nonlinear 
wave we must deduct the equilibrium angular momentum flux |a, which just 
issues from the combination of Keplerian shear and viscosity. The wave flux, we 
find, always points outwards and never exceeds that of the background viscous 
stress. The nonlinear waves never transport angular momentum inward unlike 
certain vigorously forced density waves (Borderies, Goldreich and Tremaine 
1983). 



4-4 Asymptotic description in the limit of long wavelength 



Consider Eqs (22)-(26) and suppose that their solution possesses a particularly 
long wavelength, meaning < ^ 1. In addition, let the waves oscillate near 
the epicyclic frequency so that uo = l+Odu,"^). Because we must have a = 0{1), 
it follows from Eqs (22) and (23) that u and v are 0{n^^) and that the collec- 
tive effects of pressure and viscosity (and self-gravity) are subdominant in the 
force balances. Consequently, the wave profile to leading order is determined 
by rotation and shear, a simplification that permits the momentum equation 
to be solved analytically: 

u'^ q cos'&, u'y fjL^^ q sin -& (31) 
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where the new phase variable proceeds from the transcendental equation 

e - 9o = 'd + qsin^. (32) 

Here q and are constants of the integration. The scaled amplitude q must be 
positive and less than 1 (or else neighbouring streamlines cross) and quantifies 
the rate of strain associated with the nonlinear wave, and hence the quantity 
of energy dissipated. It can be determined from Eq. (29) or at higher order in 
the asymptotic expansion. 

Note that the solution (31) is self-similar with respect to wavelength; it follows 
that our disk model can support nonlinear wavetrains of arbitrary amplitude 
and wavelength, at least formally. In the limit of large amplitude, the waves 
carry only a finite angular momentum flux, as is clear when the waveprofiles 
of (31) are substituted into (30). In particular, the coefficient of the iJ,~^ term, 
which is associated with the Reynolds stress, integrates to zero. 

The next order in the expansion appears intractable. But, if we shift to a 
pseudo-Lagrangian description we can simultaneously generalise and simplify 
the problem. At leading order collective effects arc small and the fluid el- 
ements almost exactly trace out epicycles, which may then be removed by 
averaging. We fully work out the pseudo-Lagrangian theory in Appendix A 
and summarise the main results here: the nonlinear dispersion relation and 
the equation for q. The nonlinear dispersion relation is 

u; = l + /x2Fi(g) + C>(/i=^), (33) 

in which pressure enters through the Fi function 

FM = Q-'ia - qT'^' - (34) 

Steady waves require u to be real in in Eq. (94), which returns a transcendental 
equation for q. One obtains 

{a, + fa)^F^p'\q) - 3a/?PS!j^(g) + AaqFfl^{q) = (35) 

with q — (1 — q^)~^/^ and where P^"*^ is the associated Legendre function of 
the first kind and type 3 (see Appendix A). For non- integer /3, Equation (35) 
must be solved numerically. 

4-5 Numerical computation of wavetrain solutions 

In this section solutions to (22) — (26) are computed numerically two ways: by 
the shooting and the pseudospectral methods. Our shooting algorithm inte- 
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grates the equations in 6 over one period with a 4th/5th order Runge-Kutta 
scheme with adaptive stepsizing. The algorithm converges onto the correct 
solution via the application of the four-dimensional Newton's method, which 
enforces the four boundary conditions at ^ = 27r. On the other hand, the pseu- 
dospectral method partitions the 6 domain into N points and evaluates the 
functions' spatial derivatives via excursions into Fourier space. Consequently, 
the governing equations boil down to 4A^ nonlinear algebraic equations for 
the 4 functions evaluated at the N points, plus the eigenvalue uj. Because 
of translational symmetry we may specify u — Q &X, 6 — Q which means we 
have 4A'" equations and 4A'" unknowns. The equations are solved using the 
multidimensional Newton's method. 

The input parameters are a, at, and j3 which arise from the viscosity pre- 
scription (5), and /x, the wavenumber of the wavetrain we hope to compute. 
We present results for one set of parameters only, those corresponding to an 
optical depth 1.5 (cf. Table 1), as there is little variation in the qualitative 
features as we vary the viscous parameters. Therefore, a — 0.342, ai, — 0.681, 
and 13= 1.19. 

For these fixed parameters we examined different wavenumbers We found 
travelling wave solutions for all < yUc, where jic corresponds to the marginal 
k of linear viscous overstability in Eq. (13). For the above parameters /ic = 
0.7295. If the disk was not viscously overstable (for instance, if we took the 
parameters associated with r = 0.5) then no nonlinear wavetrains could be 
found. This, of course, confirms that the wavetrains rely on the viscous oversta- 
bility to sustain them against conventional viscous loss. For ai, and a constant, 
the branch of solutions is supercritical with respect to (5. 

Wavetrains with near [x^ possess small ampUtudes, as is shown in Fig. 1. In 
this small amplitude 'linear regime' the solutions we compute connect, natu- 
rally, to the marginal linear overstable modes, which are necessarily sinusoidal. 

As ji decreases (A increases) the nonlinear amplitudes increase. For sufficiently 
long wavelengths the self-similar asymptotic profiles of (31) — (35) provide a 
good approximation to the solutions, as we can see in Fig. 2. The profiles of 
these long waves resemble those of spiral density waves, as computed by Shu 
et al. (1985a, 1985b) and Borderies et al. (1983, 1986), and also axisymmetric 
density waves in accretion disks (Fromang and Papaloizou 2007). The charac- 
teristic 'cusp '-like profile seems a recurrent feature of density waves dominated 
by the inertial forces. 

In addition to these typical profiles we show how the scaled velocity amplitude 
q varies with wavelength A = 27r//x in Fig. 3a. Recall that the actual velocity 
amplitude is equal to q/ ii (cf. Eq. (31)). By about A ~ 30 if the amplitude 
q begins to approach the value computed by the asymptotics, Eq. (35). For 
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Fig. 1. Two periods of a small amplitude nonlinear wavetrain as calculated by the 
numerical eigenvalue problem, (22)-(26). The waveprofile is plotted in the spatial co- 
ordinate of a frame comoving with the wave. The parameters correspond to r = 1.5 
from Table 1: a = 0.342, = 0.681, and f3 = 1.19. The wavenumber chosen is 
/X = 0.726 (A = 8.65) which is just below the critical ^Uc. The wavetrain profile coin- 
cides with the marginally stable overstable linear mode of the homogeneous state 
of Kepler ian shear. 




Fig. 2. Two periods of a large amplitude wavetrain; both asymptotic and numer- 
ical solution are plotted. The solid curve is the numerical solution to (23)-(26), 
the points represents the asymptotic solution (31)-(35). Parameters correspond to 
r = 1.5, as in Fig. 1, but now with fi = 0.11, so that A = 57.1. 

the parameters in which we are interested, the asymptotic result becomes an 
acceptable approximation for wavelengths above 20 — 30 H. 

In Fig. 3b we plot the angular momentum flux density over one period in the 
wavetrain as a function of A. We use the expression (30). The small ampUtude 
waves transfer negligible angular momentum, so when n Hc the angular 
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Fig. 3. Panel (a) displays the scaled velocity amplitude of the nonlinear waves q 
as a function of A = 2tt/^. Near the critical wavelength, q is almost zero and we 
are in the linear regime. From there it increases steeply and approaches a constant 
value, which agrees with the asymptotic q computed from (35), go = 0.63339. Panel 
(b) displays the total angular momentum flux. At amplitudes near criticality, the 
waves carry negligible amounts of angular momentum and the background shear 
flow dominates the transport. In contrast, larger amplitude wavetrains can carry 
appreciable amounts of up to a third of the background 3a/2. Recall a = 0.342. 

1.3 
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Fig. 4. The nonlinear wavetrain frequency u is graphed as a function of A. The 
solid curves represent the numerical solution and the points the prediction of the 
asymptotic theory, (33). 

momentum flux is approximately that of the unperturbed disk, specifically 
(3/2)a. Longer wavelength, larger amplitude wavetrains, however, can carry 
appreciable quantities above that of the background shear, but they may never 
exceed it. The angular momentum flux asymptotes to a constant value for large 
wavelengths, in agreement with the asymptotic theory. 

Lastly, in Fig. 4 the eigenvalue u is plotted against A for the same set of pa- 
rameters. The longer the wavelength the closer the frequency is to the epicyclic 
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with pressure ensuring the osciUations are a httle faster. These behaviours are 
predicted by the asymptotic nonhnear dispersion relation which appears to 
provide an adequate approximation by about A ~ 30 — 40 i?. 



4-6 Discussion of neglected physics 

The wave solutions have been computed by a relatively simple fluid model, but 
we feel that their existence and gross features are not captive to its assump- 
tions. The asymptotic analysis shows that self-gravity (a collective effect) will 
not alter the leading order wave profiles. Self-gravity will, however, alter the 
nonlinear dispersion relation (33), adding a term proportional to Also, by 
exacerbating the hnear viscous overstability it will 'shift' the curve in Fig. 3a 
to the left, which means that for a given A a wavetrain will be 'more nonlin- 
ear' (its amplitude will be greater) when self-gravity is present than when it 
is not. On the other hand, neglected thermal effects help stabilise the linear 
instability (Spahn et al. 2000) and so the curve in Fig. 3a may shift back to 
the right again. Thermal effects will also alter the pressure contribution to the 
nonlinear dispersion relation embodied by the function Fi. 

It is more difficult to ascertain the effect of vertical structure and vertical 
motion on the solutions. The linear modes may be stabilised on intermedi- 
ate scales due to their development of vertical shear (see Latter and Ogilvie 
2006b), thus shifting Fig. 3a a little to the right. But because planetary rings 
are so very thin this is probably not important. There will be additional effects 
tied to the vertical motion, the 'breathing', which allied with excluded volume 
effects at wave peaks, could lead to the vertical ejection of particles and the 
relaxation of the pressure and stress in real planetary rings ('splashing'). But 
how this plays out on moving wavecrests is difficult to forecast. The role of 
a nonlinear stress/strain relationship is also difficult to judge, though at the 
very least it will probably lead to greater effective (3s and hence shift the q 
curve to the right. 

In summary, we have demonstrated the existence and elaborated on the prop- 
erties of a family of exact nonlinear solutions to the governing equations (14)- 
(16). The solutions take the form of steady, travelling density waves propa- 
gating in either radial direction. They may also be regarded as fixed points or 
periodic orbits in the state space of the system. The members of this family 
can be distinguished by their wavenumber /i (or wavelength A) upon which 
only one restriction holds, < yUc, where fi^ is the critical wavenumber of 
linear viscous overstability. As such, an infinite number of solutions formally 
exist on arbitrarily long lengthscales; however above a certain lengthscale the 
model assumptions will certainly break down. 



19 



5 Linear stability of the nonlinear wavetrains 

The question now is: what role do the nonhnear solutions play in the evolution 
of an overstable disk? To make a start on this we need to determine the linear 
stability of the nonlinear solutions themselves. If the system is to settle on a 
wavetrain then it must be at least linearly stable. On the other hand, if it 
is linearly unstable then we are assured that the system will not settle there 
(though this does not mean it plays no role in the dynamics). We find, in 
fact, that there exists a critical wavelength Agt above which linear stability 
is assured. Wavetrain solutions with A < Agt are unstable and solutions with 
A > Ast are stable. 

These results are established by perturbing the wavetrain solution by a small 
disturbance and subsequently solving the eigenvalue problem that results for 
its growth rate. The procedure requires the solution of a Floquet boundary 
value problem which we attack analytically in the asymptotic limit of < 
/X <^ 1 (Appendix A) and numerically for general /x. 

5.1 The linear eigenvalue problem 

Consider the wavetrain associated with /i and parameters a, a^, and /3. We 
denote this basic state by Zo{9) = {(Jo,Uxo,Uyo){9), with the dependence on /i 
assumed. On this solution we superimpose a small disturbance so that 

Z^Zoi0) + Z,it,0), 

where \\Zi\\ <^ 1. Upon substituting this ansatz into the governing equations 
(l)-(2), employing the definition (19), and linearising, the system returns the 
boundary value problem: 

^tZ^ = Co{0) ^1 (36) 

where jCq is a second-order differential operator in 9 with 27r-periodic coeffi- 
cients. 

Next we assume Zi{t, 9) — e** Zi{9) which reduces (36) to a Floquet eigenvalue 
problem for the growth rate s, owing to the periodicity of Cq. In order to 
compute the spectrum we consequently make the Floquet ansatz, 

Zi = e'(^/'''^Zi(^), (37) 

where Zi is a 27r periodic function of 9, the Floquet exponent is (k/n), and k 
is the (real) wavenumber of the envelope (both s and k are different to those 
appearing in Section 3). Once substituted back into (36), one can numerically 
solve the ensuing eigenproblem for Zi and the eigenvalue s if the parameter 
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k is stipulated. 



We now present the full set of linear equations governing the initial evolution 
of the small disturbances. Instead of rendering the problem in the form of 
Eq. (36), we express it as five simpler first-order differential equations for 
Ui — Cp — Uxi, Vi = Uyi, fi, and hi. This form facilitates the numerical 
calculation. The hats will now be dropped. The five ODEs are 
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(39) 
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and these hold on the domain e G [0, 27r] subject to the periodic boundary 
conditions, (Ti(0) = mi(0) = Mi(27r), etc. Before we undertake the cal- 

culation, we need to specify the material parameters a, a^,^ /3, the wavenumber 
of the underlying wavetrain // (and hence the basic state .Zq), and lastly the 
wavenumber of the disturbance k. 



5.2 Numerical results 



The equations for the equilibrium (22)- (26) and the linear stability (38)-(42) 
are solved simultaneously using the shooting method. We restrict the param- 
eter values to those corresponding to optical depths of 1-2 and similar param- 
eters. For these values we witness little qualitative change in the results. For 
more remote values, especially ab/a small or zero, we find different stability 
behaviour but this will not be presented here. 

Once a, ab, and /? are fixed, we investigate how the various linear modes de- 
pend on the wavenumber of the underlying wavetrain /x, as a function of k. 
The reader should take care to distinguish between /i — 27r/X, the wavenumber 
of the perturbed equilibrium, and k, the envelope wavenumber of the pertur- 
bation. Normally, one need only let k vary between and /i/2; values of k 
outside this range merely reproduce modes already computed. However, in 
this section we let k take values between and This choice lets us display 
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and interpret the results in a clearer (and more familiar) way, while permitting 
us to connect the results directly with the homogeneous disk case. 

The numerical solution establishes the existence of 

• a 'generalised' viscous instability mode, which always decays at a rate pro- 
portional to ak'^; 

• two 'generalised' viscous overstability modes which can either grow or de- 
cay at a rate proportional to ak'^ and which oscillate near the epicyclic 
frequency; 

• modes which decay rapidly, at a rate of order a and which also oscillate at 
the epicyclic frequency. 

The two modulational overstable modes appear as right-going and left-going 
waves. They generally possess different growth rates because the equilibrium 
background (a nonlinear wavetrain propagating either left or right) no longer 
supports the left / right symmetry exhibited by the homogeneous disk. We are 
primarily interested in these modes, as they are the only ones which are capable 
of growing. We refer to them as 'modulational' because they exist as long 
modulations of the phase and amplitude of the equilibrium wavetrains. Thus 
the fastest growing overstable modes possess wavelengths significantly longer 
than the wavelength of the equilibrium wavetrains. 

First, in order to check the numerical results, we examined the stability of 
wavetrains of very small amplitude, i.e. those which possess a wavenumber /i 
near the critical value Hc- Recall that fic is the wavenumber of marginal vis- 
cous overstability (see Eq. (13)) and simultaneously of the very existence of 
nonlinear wavetrain solutions (sec Fig. 3). In this limit the nonlinear wavetrain 
solution, being of small amplitude, can be simply absorbed into the perturba- 
tion field and the equilibrium treated as the homogeneous state. In this limit 
the growth rates for both the 'generalised' viscous instability and 'generalised' 
viscous overstability coincide exactly with their homogeneous disk counter- 
parts (Eq. (10)), which provides the check on the numerics. The good agree- 
ment obtained between the two approaches is evident in Fig. 5 which plots 
the numerical growth rate of the modulational overstability, for the small am- 
plitude wavetrain A = 8.7, alongside the analytic growth rate obtained from 
the homogeneous problem (10). 

Next we gradually increase the wavelength of the background wave A and 
check its stabihty at each step. Some of these results we present in Figs 5 and 
6. These figures show the growth rates of the two modulational overstability 
modes as a function of k for several representative wavetrains. The model 
parameters correspond to r = 1.5 in Table 1. In the first figure, the stability 
of A = 8.7, 13.3 and 21.6 wavetrains are plotted, with the growth rate of the 
overstability mode propagating with the underlying wavetrain represented by 
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Fig. 5. Real part of the scaled growth rate s//c^ as a function of the wavenum- 
ber k for modulational overstable modes. The solid lines correspond to the modes 
propagating against the underlying wavetrain; dashed lines correspond to modes 
propagating with. Three background equilibria are selected: nonlinear wavetrains 
with wavelengths 8.7, 13.3, 21.6. The material parameters of the disk are taken from 
Table 1 for r = 1.5. In addition, we plot the stability curve for the homogeneous 
equilibrium (from Eq. (10)) as points. The amplitude of the A = 8.7 nonlinear wave- 
train is very small (being close to the critical wavelength) and thus the disk is in the 
linear regime of the homogeneous disk. Therefore the stability curve in these two 
cases agree. Larger A and wavetrain amplitudes impede modulational overstability. 
But modes propagating against the background wavetrain are more unstable than 
those propagating in the same direction. In the homogeneous disk the growth rates 
are the same. 

the dashed line, and the growth rate of the mode propagating against by the 
solid line. As is plain, increasing the A of the background (a) impedes both 
modes' growth and (b) breaks the left/right symmetry of the problem so that 
the dispersion relations of the two modes no longer coincide. In particular, 
growth of the countermoving mode is exacerbated on an intermediate range 
of k. In the limit of very small k the growth rates of the two modes are the 
same. 

From Fig. 6 we see that both modes stabilise when the underlying wavetrain is 
sufficiently long. Though the critical A at which this happens is different for the 
two modes (because of the exacerbated growth of the countermoving mode at 
intermediate k). For the parameters chosen, the modulational overstable mode 
travelling with the background stabilises when A ~ 27.0 (or 119 m, using the 
dimensions of Table 1), while the countermoving mode stabilises when A ~ 
52.8 (or 234 m). We conclude that the nonlinear wavetrains' stability criterion 
is controlled by the countermoving overstability mode. The critical wavelength 
above which a wavetrain is stable we denote by Ast- And so sufficiently long 
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Fig. 6. Continuation of Fig. 5. Real part of the scaled growth rate s/fc^ as a function 
of k for the modulational overstable modes. We present results for three background 
wavetrains: A = 39.0, A = 77.6, and the very long A = 6280 case. The material 
parameters and notation are the same as the previous figure and the range of k is 
to fi. Stars indicate the point k = ^ in both cases. The A = 39.0 graphs show the 
persistence of instability on intermediate k for the counter-propagating mode, and 
stability for its oppositely moving twin. By A = 77.6 both modes have stabilised. In 
addition, we plot the prediction of the large A asymptotic theory (see Appendix A) 
as points, and find good agreement between it and the A = 6280 case. 

wavetrains are linearly stable to all disturbances, which is in agreement with 
the asymptotic theory of Appendix A. 

Lastly, the stability of very long wavetrains was checked so as to compare with 
the asymptotic linear stability and provide a second check on the numerics. 
In Sections A. 4.1 a dispersion relation was computed giving the growth rate 
of the dominant modulational overstability mode. It is plotted with points in 
Fig. 6 against the prediction of the numerical eigensolution for a wavetrain 
with fi = 0.001. Though the agreement is good in this extreme case, the 
asymptotics provide an adequate approximation only for very long A. 



5. 3 Summary 

The linear theory is rather subtle and may repay extended study. For now, 
the essential point to take away is that there exist both linearly stable and 
linearly unstable wavetrain solutions. Wavetrains possessing a sufficiently long 
wavelength, namely A > Ast, are stable, while those that are shorter are un- 
stable. The mode which governs the stability properties is the 'modulational 
viscous overstability', in particular, the mode which travels in the opposite 
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direction to the underlying wavetrain. Those wavetrains that are unstable can 
be characterised as 'saddle points' in the phase space because they possess the 
unstable eigenmodes (the modulational overstabilities, which grow slowly, no 
faster than ak'^), and a set of stable eigenmodes (some of which decay speedily, 
as a). 

This general behaviour is repeated for all realistic parameters we tried. If we 
turn to the parameters listed in Table 1, a ring with r = 1 yields Agt = 50.2 
(165 m), a ring with r = 1.5 yields Agt = 52.8 (234 m), and a ring with 
T — 2 yields Agt = 117 (640 m). The much longer value in the last case is 
due to the system's sensitivity to the f3 parameter. For fixed a and 05 we find 
that Ast is a steeply increasing function of f3, a trend that is summarised in 
Table 3. The critical lengths at r = 1-1.5 compare suggestively to the scales 
of quasi-periodic microstructure observed by Cassini, which are some 150-220 
m (Colwell et al. 2007, Thomson et al. 2007). But it should be noted that 
important physics remains missing and the stability estimates are subject to 
some revision. In particular, the results for r = 2, and greater f3 generally, are 
probably blemished by neglect of certain physical processes (self-gravity, for 
example). Nevertheless, the overall consistency is encouraging. 





Ast (in H) 


1.1 


41.9 


1.19 


52.8 


1.3 


71.0 


1.4 


91.7 


1.5 


121 



Table 3 

The critical A of stability for nonlinear wavetrains as a function of (3. The other 
parameters a and 05 are held fixed: a = 0.342 and aj, = 0.681. The case (3 = 1.19 
corresponds to an optical depth of 1.5, larger values of /3 correspond to more optically 
thick disks. For instance, the (3 = 1.5 case is associated with r = 2 (Schmidt et 
al. 2001). 



6 General nonlinear dynamics 

Armed with the linear stability results, we return to the question posed earlier: 
what role does this family of 'fixed points' play in the general dynamics of an 
ovcrstable disk? The simplest outcome would be for an overstable disk to 
migrate from the homogeneous state to the shortest stable wavetrain solution 
(with A = Agt); but this is far from guaranteed. The dynamics is likely to be 
more complicated, with the system flirting with a number of stable nonlinear 
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solutions and thus exhibiting irregular time and space dependent variations. In 
this section we explore this behaviour with some schematic ideas and analogies. 
We will not bury ourselves in technical details; this we leave for later, when 
we have a more complete set of numerical results at hand. 

First, we offer an explanation of the cascade of power to longer lengths that 
has been observed in nonlinear hydrodynamical and A^-body simulations. The 
analysis is framed in the language of dynamical systems theory and though 
'hand-wavey' should be straightforward to check numerically. Next we exam- 
ine the nonlinear interactions between linearly stable wavetrains of different 
wavelength, by computing slow variations in the wavetrains' phase. It can 
be shown that the disk supports modulational 'shock' and 'source' structures, 
whereby the wavenumber and/or amplitude of a wavetrain undergoes a discon- 
tinuity. Lastly, we draw analogies between a viscously overstable disk and the 
dynamics of the complex Ginzburg-Landau equation, and subsequently make 
a few additional predictions about the likely nonlinear behaviours associated 
with overstable disks. 



6.1 The cascade to longer lengthscales 

In Fig. 7 we have drawn a schematic diagram of the state space of the over- 
stable disk. The state space is infinite-dimensional, and, needless to say, some- 
what difficult to visualise, but our crude two-dimensional projection raises a 
few interesting points and predictions. This representation of the dynamics 
treats the time-evolution of the state vector Z as a trajectory, tracing the 
smooth transition of the system from one state to another. Equilibria appear 
as invariant fixed points: trajectories that begin there stay there, and their 
linear stability decides if trajectories that pass infinitesimally close fall into 
the fixed point or are repelled. The state of homogeneous Keplerian shear is 
such a fixed point, and when it is overstable it behaves like a saddle. Nearby 
trajectories will be defiected away from the point along those directions as- 
sociated with the overstability modes. Simultaneously, the same trajectories 
will be drawn towards the fixed point along those directions associated with 
the stable modes. 

A nonlinear wavetrain solution should appear as a periodic orbit, but for 
clarity we represent it as a fixed point in Fig. 7; our argument is unchanged. 
The entire family of these nonlinear wave solutions plots out a semi-infinite, 
one-dimensional curve in the state space, which we parametrise by A, the 
wavelength. The foot of this curve is the homogeneous state. In Fig. 7 the 
thick curve represents this continuous family of solutions, and the thinner 
arrowed curves represent possible trajectories of the system. 
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Homogeneous Disk 



Fig. 7. A schematic drawing showing a two-dimensional projection of the overstable 
disk's state space, and the upward cascade to longer lengthscales. The thick curve 
represents the continuous family of nonlinear wavetrains, parametrised by A, and the 
thin curves represent typical trajectories. Also plotted are the state of homogeneous 
Keplerian shear and the first stable nonlinear wavetrain, for which A = Ast- The 
basic idea is that trajectories repeatedly pass from the overstable manifold of one 
fixed point to the stable manifold of a neighbouring longer wavelength fixed point. 
This process ends, or at least undergoes a qualitative change, once the trajectory 
enters the vicinity of the first stable fixed point. 

Waves are unstable if A < Ast, where Ast is some critical wavelength, while 
longer waves are stable, and like the overstable homogeneous state, the unsta- 
ble wave solutions are saddles, each possessing a slow overstable manifold and 
a fast stable manifold. Trajectories near the unstable segment of the curve of 
fixed points will be drawn towards the curve along the nearest stable mani- 
fold; at the same time trajectories are directed 'up' the curve by the overstable 
manifold. Because the time-scale of the attraction is faster than that of repul- 
sion, trajectories will be bunched, or 'focused', about the curve of fixed points 
even as they travel up it. Close to the branch of invariant fixed points this 
migration may be facilitated by the the linear 'translation mode' with s = 
and eigenvalue dgZ^ (see next subsection). Thus the thick curve in Fig. 7 is 
a one-dimensional 'spine' around which collects a nest of state-space trajecto- 
ries. Such behaviour is typical of dissipative systems, in which the long time 
dynamics is governed by a lower dimensional subset of the state space (e.g. 
Robinson 2001). 

The upward drift will end, or at least change qualitatively, once a trajectory 
reaches the vicinity of the first stable fixed point. Because larger A fixed points 
are, necessarily, associated with longer wavelengths this upward migration 
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should coincide with the movement of power to longer and longer scales. 



It is not immediately obvious what happens once a trajectory nears the first 
stable fixed point, nor how a trajectory behaves if it begins near the stable 
portion of the solution branch. The simplest outcome is for the system to settle 
on the first stable fixed point available (Agt for example), but this is far from 
assured. The basin of attraction of a stable fixed point may be very small, and 
we expect this to be the case for most parameters. It is more likely that the 
system will 'wander around' the set of stable solutions, yet never fall upon any 
particular one. Physically this 'wandering' will appear as temporal and spatial 
modulations upon a bed of nonlinear waves, with the modulations manifesting 
on lengthscales longer than that of the underlying wavetrains. It is not clear 
which set of A the system will select for the underlying waves, or whether this 
set gradually changes (perhaps growing larger and larger with time). Moreover, 
the situation may be complicated by a chaotic attractor (as can be the case in 
the complex Ginzburg-Landau equation). That said, it is plausible that small 
initial conditions will yield a band of saturated wavelengths near Agt. 



These conjectures are straightforward to check with numerical simulations, 
and we are currently undertaking this work. The only published simulation, 
that of Schmit and Tscharnuter (1999), exhibited an upward cascade which 
in the nonself-gravitating case had not yet halted after 10, 000 orbits. By that 
stage most power was situated on a scale of some lOOH. According to the 
linear stability calculations of the previous section, using the same parameters 
as Schmit and Tscharnuter, we find Agt ~ 60H; so the simulation definitely 
passed the critical wavelength of the hnear analysis. Recall, however, that these 
simulations do not support global nonlinear wavetrain solutions on account of 
the boundary conditions and this negates much of the interpretation of this 
section. In their simulations, trains of travelling waves can emerge locally but 
given sufficient time (and 10,000 orbits is more than enough) these pulses will 
traverse the domain, encounter the boundaries, and reflect backwards and 
interfere with themselves. Such nonlinear interactions are not captured in our 
model and perhaps induce an injection of power to longer scales than what 
we would expect. 



Finally, one should note that in our interpretation the halting of the upward 
cascade does not require self-gravity, which is an effect Schmit and Tscharnuter 
emphasise. In their simulations its inclusion halts the upward cascade. But it 
is probable that the physical mechanism is different to that we have sketched. 
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6.2 Modulated wavetrains: weak shocks, sources, and cellular dynamics 



This and the next section ascertains the basic nonhnear dynamics of wavetrain 
modulations. This will give us some handle on the behaviour of the system 
once it nears the set of stable fixed points described in Fig. 7. The problem is 
a difficult one, and not generally amenable to analytic techniques. That said, 
a useful result can be derived if it is assumed that the modulations in question 
vary slowly in time and space. When this is the case, the phase perturbations 
are governed by the Burgers equation, which in turn suggests that the radial 
domain may fracture into regions ('cells') of different wavenumber bounded 
by two kinds of interface: 'weak shocks' and 'weak sources', where the phase 
undergoes rapid change. The actual dynamics is perhaps even more compli- 
cated than this ideal picture but it remains a useful paradigm to understand 
the larger-scale irregular variations observed. It is also worth remarking that 
the theory is akin to the attractive idea proposed by Tremaine (Araki and 
Tremainc 1986, Tremaine 2003) whereby disk structure corresponds to 'jams' 
— the key difference is that in the Tremaine model the jams correspond to dis- 
continuities in shear, and in our model the jams correspond to discontinuities 
in the wavenumber of background waves. 

We now very briefly describe the derivation of the Burgers equation. The proof 
is lengthy and tedious so more details can be found in Appendix B. It is essen- 
tially a multiple-scales analysis which applies the methodology of Howard and 
Kopell (1977) and Doelman et al. (2009) who examined modulated wavetrains 
in reaction-diffusion systems. 



6.2.1 The Burgers equation 

Consider a linearly stable wavetrain of fi and uj = uj{^) when ^ is small. 
We represent it by the vector Zq{6 ; /x), where we have made the dependence 
of Zq on wavenumber /i explicit in contrast to Section 4. The wavetrain is 
assumed stable; therefore, associated with it are a set of decaying linear modes 
possessing negative growth rates. Of these we select the modulational viscous 
instability mode and denote its growth rate by s{k) where k is its wavenumber. 
From Section 5 and Eq. (103), this mode always decays as k"^. 

Slowly varying modulations of Zq are sought. First we establish the long length 
and time scales characteristic of this modulation. For a small dimensionless 
parameter < 5 ^ 1, define new (slow) variables 

X^5{x-Cgt), T^S^t, 
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where Cg is the group velocity of the underlying wavetrain: 

_ duj 

Note that we have chosen a spatial frame moving with the group velocity. 
Wavetrain modulations travel at c^, not the phase velocity Cp. 

Next a perturbation of the wavetrains' phase is introduced, Q{X,T). The 
associated variation in wavenumber is 9x©. We now construct a solution of 
the form 

Z = Z^{e + Q- ^ + 5Qx) + 6^ Zi{e,X,T), (43) 

where a subscript X (or T) indicates partial differentiation. Our strategy is to 
derive an equation for the phase modulation © so that the ansatz (43) satisfies 
the nonlinear equations to the leading orders in 5. 

The ansatz (43) is substituted into (14)-(18) and these equations are expanded 
in powers of the small parameter 5. At leading order 5° the balance becomes 
the nonlinear eigenproblem for the underlying wavetrain profile, Eqs (23)-(26). 
At the next order 5^ the balance is a simple identity: the first ji derivative of 
the leading order equations. At order 5'^ we obtain 

Co Zi = -{deZo) Qt + A{e) Qxx - B{e) 9^, (44) 

where Cq is the linear operator introduced in Section 5, the A and B are 

vectors that depend only on the underlying wavetrain and so depend only 
on and 9. Their expressions are complicated and we omit them here (see 
Appendix B). The order 5^ equation need only yield a solvability condition: the 
right side of Eq. (44) must he in the range of Cq (the Predholm alternative) . 
This can be assured if the inner product of the right side with the adjoint 
solution of Cq is zero. After some laborious manipulation, the condition is 
equivalent to the Burgers equation, 

Kt + KKx + |s"(0) Kxx = 0, (45) 

where K = Qx is the wavenumber modulation, and 

The 'phase diffusion coefficient' is — s"(0)/2 and is associated with the de- 
caying modulational viscous instability. Localised disturbances which 'bunch 

up' the crests of wavetrains will relax according to this mode. The 'advective 
coefficient' is uj" and thus equal to the group velocity's rate of change with /i. 
It controls the steepening of fronts or shocks in K. If Cg is independent of \i 
there is no wave steepening. How one derives s"(0) from A and a;"(A;) from B 
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is nontrivial and is left to Appendix B. There we also show how to compute 
the coefficients in terms of the material properties of the disk, a, 0:5, and (5. 

6.2.2 Weak shocks and sources 

For our purposes, equation (45) predicts two basic behaviours. Localised solu- 
tions K{X, T) will decay to zero as T ^ 00 while propagating at a speed equal 
to the group velocity (Whitham, 1974, Doelman et al. 2009). This behaviour 
simply expresses the linear stability of the underlying wavetrain. Nonlocalised 
solutions, on the other hand, can manifest as viscous Lax shocks. For instance, 
suppose that we require 

K — > K+ as X — > 00, 
K — > K_ as X — > —00, 

where and K_ are two real constants. Then Eq. (45) admits a travelling 
front or shock solution i^(X, T) = K^h^X — CshT) when 

u;"{i^)(K+- K_) <0. (46) 

From Fig. 4 or from Eqs (33) and (34) we have uj" > and thus < K^. 
In addition, the speed of the front is given by 

Csh=y(/^) {K+ + K_), 

which is the Rankine-Hugoniot condition in the moving frame. The character- 
istics associated with the solution point towards the front interface, as do the 
motions of the individual wavecrests on its each side. Solutions for which the 
characteristics and wavecrests point away do not satisfy (46) and correspond 
to rarefaction waves or 'sources'. Fig. 8 presents a cartoon of a nonlinear wave- 
train exhibiting both a shock and a source at which points the wavenumber 
(and amplitude) rapidly change. 

6.2.3 Cellular structure 

In general reaction-diffusion systems, widely-spaced sources and shocks can 
partition the domain into distinct regions inhabited by wavetrains of different 
H (Aranson and Kramer 2002, Doelman et al. 2009). Figure 8 roughly sketches 
a portion of a domain so decomposed. Similar phenomena should emerge in 
overstable disks, and it is possible that observed radial structure on interme- 
diate scales in Saturn's rings corresponds to the cellular patterns constructed 
from such wave defects. 

Shock and souce interfaces will possess their own slow dynamics which compels 



31 




source shock 

radius 



Fig. 8. A simple cartoon representing a wavetrain with a shock and a source solution. 
The arrows represent the direction that the individual wavecrests propagate (the 
Cp), while K represents the modulational wavenumber. Here Ki < K3 < K2 are 
constants, and thus the defects 'glue' together patches of different wavenumber 
(and hence amplitude). On sufficiently large scales a wavetrain may decompose into 
a number of subtrains each possessing varying wavenumber and each enclosed by a 
shock and a source. These interfaces will possess their own slow dynamics — they 
may wander around the domain, annihilate one another, generate new defects, etc. 

them to drift relative to each other and undergo various kinds of interactions. 
The former motions may be represented by a simple one-dimensional 'cellular' 
dynamical system that evolves the interfaces according to a set of ordinary 
differential equations (Elphick et al. 1990, Ei 2002, Doelman et al. 2009). Each 
interface is treated as a particle that weakly interacts with its neighbours via a 
'potential' and thus the model bears a similarity to a one-dimensional A^-body 
problem. 

Unfortunately, it is difficult to produce analytic estimates for the character- 
istic lengthscales of the cellular structures so formed. Thus comparison with 
Cassini data must wait. In any case, the model assumes only weak inter- 
actions between the interfaces, but full numerical simulations of the com- 
plex Ginzburg-Landau equation typically exhibit more complicated behaviour 
(Popp et al. 1993, Chate 1994). We briefly survey these phenomena in the 
next subsection. 

6.2.4 The complex Ginzburg-Landau equation 

The complex Ginzburg-Landau equation in an extended domain boasts a sub- 
stantial literature, equal to the diverse behaviour it exhibits. This body of work 
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may guide the interpretation of future nonlinear simulations of the overstabil- 
ity. The complex Ginzburg-Landau system crystallises the essential dynamics 
of general nonlinear wave phenomenon in a broad range of physical situations 
(Aranson and Kramer 2002) and which we expect nonlinear waves in planetary 
rings to share Q 



The one-dimensional version of the equation reads 

dtA = A + {l + %cx)dlA - (1 + iC2)\A\^ A, 

where A is a complex function and c\ denotes the linear dispersion and C2 
nonlinear dispersion. The function A typically represents the complex ampli- 
tude of some wave phenomenon, and is not the same as the phase modulation 
O which we studied earlier — the latter, though, can be derived from A. The 
equation gives rise to steady, nonlinear wavetrains, various wave instabilities, 
and 'defects' (such as source and shock structures) that can 'glue' together 
different wavetrains, in addition to different types of spatio-temporal chaos 
(Aranson and Kramer 2002). 



Of the various regimes exhibited, behaviour associated with the so-called 'in- 
termittency regime' is of most interest. Here the state space admits not only 
stable nonlinear wavetrain solutions but also a 'defect-turbulent attractor' 
(Chate 1994). The role of the latter is to govern the irregular dynamics of 
the shocks and sources that separate the different regions of stable travelling 
waves. A number of studies have explored the rich and chaotic dynamics of the 
interfaces (Popp et al. 1993, Chate 1994, Ipsen and Hecke 2001): for instance, 
source and shock structures may wander around the domain, annihilate each 
other, and give birth to new shock and source pairs. The long-time outcome 
is also variable, and depends closely on the parameters. For certain choices 
of Ci and C2, the final state may be 'glassy' (a static arrangement of defects), 
the defects may evaporate altogether, or the birth and annihilation process 
continues indefinitely. Another interesting regime is that of 'phase chaos' in 
which the wavenumber of the nonlinear wavetrains undergo turbulent fluctua- 
tions. All or some of this behaviour we expect to carry over to the dynamics of 
the overstable disk. At the very least, the complex Ginzburg-Landau equation 
provides a template with which to interpret simulations of overstable disks 
and the ISS observations of irregular structure in the B-ring. 



In fact, the complex Ginzburg-Landau equation has been derived from a weakly 
nonlinear analysis of the overstability in a simple hydrodynamic model (J. Schmidt, 
private communication) . 
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7 Conclusion 



This paper has investigated the dynamics of an overstable fluid disk by ex- 
ploiting the properties of its nonhnear invariant solutions. These structures 
manifest as axisymmetric uniform travelling wavetrains and form a continuous 
family of solutions parametrised by their wavelength A. In Section 4 and Ap- 
pendix A we computed the members of this family and outlined their physical 
characteristics. In particular, we flnd that inertial forces dominate their leading 
order force balance, and so the fluid undergoes approximately epicyclic mo- 
tion. The wave amplitudes, on the other hand, are determined from an energy 
balance that counterpoises linear viscous overstability and viscous dissipation. 
In Section 5 the wavetrains' stability to small perturbations is inspected. It 
turns out that only solutions with wavelengths above a critical value Agt are 
linearly stable. If the viscous parameters take values associated with Saturn's 
B-ring, Ast ~ 50 — 100 H. Finally, in Section 6 we speculate on the role of 
the nonlinear wavetrains in the general dynamics, arguing that the system 
will linger 'near' the branch of linearly stable invariant solutions. Physically, 
this phase-space behaviour will translate to a bed of nonlinear waves suffering 
large-scale modulations or defects, the latter partitioning the domain into a 
cellular grid of wavetrains of different A. 

Our basic hypothesis is that quasi-periodic microstructure on scales of order 
100m revealed by the UVIS and RSS correspond to the underlying bed of 
nonlinear waves. Meanwhile, irregular structure on larger scales, of 500m to 
10km may be a manifestation of the wavetrains' spatio-temporal modulations, 
arising as the system flirts with a number of stable wavetrain solutions but 
is unable to settle on any. The individual waves arc beneath the ISS cam- 
eras' resolution, but larger-scale variations in their amplitude and wavelength 
should give rise to different optical properties which the cameras can register. 
In order to fully establish the latter point, however, the optical properties of 
the nonlinear waves need to be understood. The dynamical and photometric 
modelling techniques that have been developed to understand the observa- 
tional properties of gravitational wakes could be applied to this end (French 
et al. 2007, ColweU et al. 2007). 

This study has adopted some rather strong modelling assumptions and ne- 
glected some important physics. The issue requiring the most urgent attention 
is self-gravity, as it will certainly play a role in the formation and stability of 
the nonlinear waves studied. Its effects we examine in a following paper. The 
ID isothermal fluid model we have employed is simple, but we doubt that the 
basic physics it has revealed is an artefact of its simplicity. That said, on at 
least two counts it can be improved: first, its treatment of the viscous stress, 
and second, its neglect of vertical motion. The Navier-Stokes stress prescrip- 
tion is not ideal in a scenario of densely-packed and infrequently colliding 
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particles, where the coUisional stress dominates. A detailed kinetic treatment 
is perhaps the best way to captm^e the nonlinear dependence on strain associ- 
ated with a coUisional stress (Latter and Ogilvie 2008). Vertical motions, such 
as 'breathing' or 'splashing', certainly accompany both the linear modes of the 
viscous overstability and its nonhnear development (Salo et al. 2001, Latter 
and Ogilvie 2006b). The importance of this effect is difficult to judge; and may 
only be resolved by 2D eigensolutions and simulations, both formidable tasks. 

Lastly, we briefly sketch out future work. Once the role of self-gravity on the 
nonlinear waves has been established, large-scale hydrodynamical simulations 
of an overstable ring must be undertaken to test our various conclusions. These 
should also address the influence of the boundary conditions, which will help 
us connect the work here to that of Schmit and Tscharnuter (1999). On the 
other hand, the dense-gas kinetic equations may be simulated, using perhaps 
the continuum model derived by Latter and Ogilvie (2008). These then could 
provide a more direct comparison with A^-body simulations and possibly to 
the observations of finescale structure themselves. 
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8 Appendix A: Asymptotic description of ring dynamics in the 
limit of small collective effects 

The analysis presented in this appendix resembles earlier work on nonlinear 
waves in disks. Several papers by Shu and collaborators formulate equations for 
nonlinear spiral density waves in planetary rings and accretion discs (the most 
pertinent being Shu ct al. 1985a, 1985b). They include self-gravity, pressure 
and viscous (or kinetic) stresses. However, they do not consider time-evolution 
or viscous overstability because they are interested in forced rather than free 
waves. The streamline formalism of Borderies et al. (1983, 1986) is closely 
related. They derive their equations using methods of celestial mechanics but 
the outcome is similar. Because they consider discrete streamlines rather than 
a continuous medium, they obtain coupled ODEs rather than PDEs which 
means that density waves are not so clearly manifested in their approach. Our 
analysis is intentionally restricted to axisymmetry and local geometry, which 
slightly clarifies matters and allows viscous overstabihty to take centre stage. 
The physics of self-gravity, nonaxisymmetry, and the gravitational infiuence 
of external moons may be included later. 

8.1 Basic equations 

We consider axisymmetric motions in a shearing sheet with no vertical exten- 
sion. The sheet is Keplerian and units are chosen such that Q — 1. The basic 
equations are 

Du-2v^ eT^, (47) 
Bv + ]-u^ eTy, (48) 
Dcr + ad^u = 0, (49) 
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where D = dt + udx is the Lagrangian derivative, and Ty are the radial and 
azimuthal accelerations due to the weak collective effects of pressure, viscosity, 
and self-gravity, and e is an ordering parameter. Note that the notation here 
is different to the main body of the paper: u denotes the radial velocity and v 
the azimuthal velocity. 



We anticipate that the dominant motion of fluid elements is an epicyclic os- 
cillation of the form 

x^i + a{i,t)cos[t-h{i,t)l (50) 

where a and h are the epicyclic amplitude and phase, and ^ is an approximate 
Lagrangian variable. It is not an exact Lagrangian variable because the exact 
motion of a fluid element in a viscous disk departs marginally from epicycles. 
In fact, we use this equation to define C,{x, t), and then determine how a and b 
must evolve so that this equation continues to describe the dominant motion 
of the fluid. This constraint is meaningful only if a and b evolve slowly in time, 
a constraint enforced by the asymptotic analysis below. We write a^ and at 
to represent da/d^ and da/dt, respectively. We also work with the complex 
epicyclic amplitude 

A = ae'\ (51) 

The relation between x and ^ is one-to-one provided that 1^4^ | < 1. Interpen- 
etration of neighbouring fluid elements (or streamlines) occurs if \A^\ > 1. 



We now change variables from {x,t) to {C,t), flnding 

Du-2v^ eT^, (52) 

Bv + ^u^ eTy, (53) 

Dcr + aJ-^d^u = 0, (54) 

where now 

D = 9t + uJ-^d^, (55) 

u — u — X, (56) 

(Ox \ 

—j =l + a^ cos{t -b) + ab^ sm{t - b) , (57) 

(dx \ 
—j = at cos{t — b) — a{l — bt) sm{t — b) . (58) 

(Note that dt now denotes a derivative at constant J is the Jacobian 
(stretching factor) of the transformation between x and ^, while x is the radial 
velocity of a particle having constant ^. If ^ were an exact Lagrangian variable, 
then u would equal x, therefore u is the radial velocity relative to the moving 
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coordinate system. Note the identity 

dtJ = d^x. 



(59) 



8.2 Asymptotic solution 

We propose that a and b evolve, due to weak collective effects, on a timescale 
that is long (by a factor of order e~^) compared to the orbital timescale. Thus 
a and b are considered as functions of where t — et is a slow time 

variable, and at becomes ear, etc. 

The solution is expanded in the form 

u — eui + • • • , 

V = vq + evi -\ , 

a = ao + eai -\ , 

where ui, etc., are functions of (^,t,r). (Note that these functions are dif- 
ferent to those appearing in Section 5.1 though the notation is the same.) 
The multiple-scale analysis is necessary in order to separate the fast epicyclic 
oscillation from the slow evolution of the epicyclic amplitude and mass distri- 
bution. All the variables should be periodic functions of t with period 27r; any 
secular evolution is described by the dependence on r. 

Note that a and b are not expanded in powers of e, because they represent 
only the leading-order motion. Thus J(^, t, r) is a quantity of order unity and 
does not require expansion. On the other hand x contains terms of orders e° 
and e^. The relation between u and u is 

u — u = X = — asin(t — b) + e[ar cos(t — 6) -|- ab-j- sin(t — b)] = xq + exi. 

The identity (59) becomes 

dtJ = d^xo, drJ = d^xi. (60) 

At leading order we have the epicyclic motion 

Xq — —asm.{t — b), Vq = —-acos{t — b). (61) 
This satisfies the equations of motion (52) and (53) at 0(1). The equation of 
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mass conservation (54) at gives 

dtao + aoJ~^d^Xo = 0. (62) 

Multiplying by J and using the identity (60), we obtain 

dt{aoJ) = 0, (63) 

which expresses conservation of mass on the epicyclic timescale. We may write 
(To J = E(^,r). 

Our task now is to deduce evolutionary equations for the epicyclic amplitude 

A{^,r) and the epicycle-averaged surface density E(^,r). (Note that, since 
dx — J at fixed E is really a density with respect to the ^ coordinate.) 



8.2.1 Evolution equations 

The equation of mass conservation at 0{e) is 

dtai + {dr + UiJ~^d^)ao + aoJ~^d^{ui + Xi) + a^J'^d^Xo = 0. (64) 
Multiplying by J and using the identity (60), we find 

dt{(TiJ) + dr{aoJ) + d^{aoUi) = 0. 
We now apply the epicyclic averaging operation 




to obtain 

+ d^iEV) = 0. (65) 

where V{^, r) is a mean radial velocity defined by the averaged mass fiux 
density 

EV = {aoui) = E{uiJ-^). (66) 

What we have done here is to obtain an evolutionary equation for E by ap- 
plying a solvability condition that eliminates the unknown a^. 

The equations of motion at 0(e) are 

dtiui + xi) + (dr + uiJ-'d^)xo - 2vi = T,, (67) 
dtvi + {dr + u^J-^d^)vQ + ^(ki + ii) = Ty. (68) 

Our aim is to eliminate the unknowns Ui and Vi by obtaining a solvability 
condition that will provide an evolutionary equation for the epicyclic ampli- 
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tude A. We will show how this can be obtained in a systematic manner. Let 
us rewrite the equations, noting that d^Vo = (1 — J)/2 and St-Vq — —Xi/2, in 
the form 

dtui + -j^tJ - 2f 1 = = 7; - dtxi - drXo, 

The latter (angular momentum) equation, when averaged, gives us immedi- 
ately 

V^2{Ty). (69) 
Now eliminate Ui between the two equations to obtain 

-2Jvi - 2dt{J^dtVi) = Jf, - 2dt{JX) 

On the left-hand side we have a linear operator in self-adjoint form acting on 
vi. It can be verified that J~^e'* is a null eigenf unction. The corresponding 
solvability condition is 

{J-^e'yf, - 2dt{JX)]) = 0. (70) 
After an integration by parts and use of the identity 

we have 

+ 2i{e'' + A^)Ty) = 0, 

and so 

drA + Vd^A = {(^'{iT, - 2Ty)). (71) 
In summary, we have derived the evolutionary equations 

+ d^{Ty) = 0, (72) 

drA + Vd^A = {e''{iT, - 2Ty)), (73) 

with 

V = 2{Ty). (74) 

The evolution of the mass distribution and the nonlinear waves are coupled 
together. We can now remove the ordering parameter and replace dr by dt 
here if desired. However, it is important to retain the conceptual distinction 
between the fast and slow timescales. 
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8.3 Evaluation of the force averages 



8.3.1 Stress tensors in general 

If the forces derive from a (two-dimensional/ vertically integrated) stress ten- 
sor, then 

Tx — 9xf, Ty dxh, 
a a 

where / and h are as given in Sections 3 and 4 in the case of a fluid with an 
isotropic pressure and a Navier-Stokes viscous stress. Self-gravity is omitted. 
To calculate the stresses to leading order, we replace a with (Jq here and use the 
leading-order velocity field where required. Also, noting that aQ^dx — T,~^d^ 
and that E is independent of the fast time t, we have 

{e^'iiTx - 2Ty)) = ~d^{e'\if - 2h)), (75) 
{Ty)^~d^{h). (76) 



8.3.2 Pressure 

We first compute the pressure contribution to the evolutionary equations by 
setting the viscous stress to zero. In the case of an isothermal gas with (two- 
dimensional/ vertically integrated) pressure p = v^a, we have 

f^v',EJ-\ h^O. (77) 

We must then evaluate (e**J~^). Let 

= (a^ + iab^)e''' = qe''^. (78) 

Here < g < 1 is the usual nonlinearity parameter for density waves in 
planetary rings (see Borderies, Goldreich, and Tremaine 1983, 1985, 1986). 
Then 

J = 1 + Rc[A^e-'^] = 1 + Re[qe-'^^~'f'^] = 1 + qcos{t - <j)) ^ 1 + qcos0 

and 

/e''j-^) = e'^{e'\l + qcos9)-^) = e"^— / cos^(l + gcos^)"^ rf^, 

ZTT Jo 

where 9 — t — (f). This is straightforward to integrate and the evolutionary 
equations with pressure alone become 

dt^ = 0, (79) 
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d,A=-^d^[F,{q)J:A^], 



(80) 



with 



Fi{q)=q-' (1-g 



This is a type of nonhnear Schrodinger equation. Note that q^ 



8.3.3 Viscosity 

Next we compute the contribution from viscosity alone, in the absence of 
pressure. Then 

and we adopt the viscosity laws 

where is a constant reference density. To leading order and in terms of the 
variables Q = ge**^ and 9 = t — 4> introduced previously, 

/ = vl {ah + |a)a, i^—j r^-^q sin 9, (84) 

h^vlaa^y—j J-^-^(| + 2gcos^). (85) 

The weighted averages of these relations, which appear in Eqs (72)-(74), are 
complicated integral functions which may be rewritten in terms of associated 
Legendre functions. In summary, 

{e^\if - 2h)) = -ay, l^-J FM (86) 
{h) = -ay FM. (87) 

where we have introduced two new F functions defined through 

^^^^^ = /j(lt' J) q {(«^ + f «)|P^'^(^) - 3a/3P2^(g) +4agPg^(g)| 

(88) 

^^^^^ ^ "TTI {^^^ + " 4gPS5^(g)} , (89) 



h = -va{-\ + ^a^-y). 



(82) 



(83) 
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where g = (1 — q^) ^Z^, a notational convenience, and P^"*-* is the associated 
Legendre function of the first kind and type 3. It is defined through 



r(m — 7]) 



TTT{r]) Jo 



cos{m9) z + V z"^ — 1 cos 9 



-,,-1 



de 



(see|http: / /functions.wolfram.com/HypergeometricFunctions/LegendreP3General ) 



8.3.4 Combined evolutionary equations 



Dimensions are adopted so that f ^ = 1 and a* = 1. Finally, we arrive at a 
problem of the form 

+ = 0, (90) 



dtA + Vd^A=-d^[[tFM + ^^F2{q)\ SA^} 



where 



(91) 



(92) 



and where the Fj's are the real functions presented in Eqs (81) and (88)-(89). 



The behaviour of the Fi controls to a large extent the evolution of the system. 
When solutions are of small amplitude and q is near 0, the three F functions 
are approximately constant because P^'"^(g) oc g™. In particular ^3(0) = —3a 
and 

^2(0) = |[3afc-2a-9a/3]. (93) 

So, we have linear viscous overstability (antidiffusion) when -F2(0) < (cf. 
Eq.(12)), and viscous diffusion when -F2(0) > 0. On the other hand, if the disk 
is overstable then F2 must change sign as q increases because when q ^ 1 
F2 is positive. This suggests that the instability saturates once the solution is 
sufficiently nonlinear. The function F3 also changes sign at a critical value of g, 
a situation which corresponds to angular momentum flux reversal. It appears, 
for all realistic parameters, however, that viscous overstability saturates at a 
q below the critical value of flux reversal. In the opposite limit of g — 1, F2 
diverges as (1 — q^)~l^~^l'^ because viscous diffusion is trying to combat the 
imminent crossing of neighbouring streamlines. 



8.4 Nonlinear travelling waves 



The governing set of equations (90)-(92) appears rather complicated, never- 
theless it admits exact plane wave solutions of the form 



A = Aoe 



i^^—iuit 
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where fj, and cu arc real so that q = qo — |/^^o| = constant. The solution also 
requires E = constant and V = 0. These represent uniform travelling waves in 
a homogeneous disk and provide an asymptotic description of the nonlinear 
waves studied in the main body of the paper. 



Steady wavetrains require 

-icu^ -fi^[iFi{qo) + F2(go)], (94) 

which is possible only if ^2(^0) = 0. Therefore the amphtude of the nonlinear 
waves (in the sense of q) corresponds to the condition of marginal overstability. 
We have go — qc, where qc is the critical q at which ^2(5) reverses sign. (If 
-^2(0) > 0, there is no linear overstability and no travelling wave.) However, the 
amplitude in the sense of Aq is not limited. The nonlinear dispersion relation 
fixing uj is consequently controlled by pressure: uj = iJ? Fi(go)- 



8.4.. 1 Linear stability 



Consider small perturbations upon the uniform travelling wave solution just 
described. The small perturbations are denoted by 5E, 5 A, etc. Their linearised 
equations are 

dt5T. + T.d^5V = 0, (95) 

dM + 6Vd,A = + d,{6q'A,) 

- + iF^d^ [^"^i) + 'F^^A^^ (96) 

<^^=||W + (/3 + im(f), (97) 

with 

Sq^ = SA^Al + A^SAl 

Solutions are of the form 

5E = j:(ciE + clE*), 
5V = C2E + c*E*, 
SA^A{c3E + clE*), 

where 

E — ex-p{st + ik^), 
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s is a complex growth rate, k is a real wavenumber and the q are complex 
coefficients. This allows us to write 



Sq' = q' 



C3 1 



C4 1 



E 



+ 



cU 1 + - + cl 1 - - 



E* 



And now we have a fourth order algebraic system for the coefficients q, 



sci = —ikc2 
C2 = ikPs 



03(1 + -) +04(1-- 



+ iA;(/3 + l)F3Ci, 



(s - ia;)c3 + i/i.C2 = -(F2 + iFi){fi + /c)[c3(/x + /c) + C4(/i - A;)] 

— iFifikci — i{n + k^FiCs, 

[s + ia;)c4 - iiiC2 = -(-F2 - iFi){l^ - k)[c-i{^, + k) + C4{fi - k)] 

— iFinkci + i{n — k^FiCi, 



(98) 
(99) 

(100) 
(101) 



where the Fi arc evaluated at go and Fj = q^dFi/dq^ evaluated at go- A 
complicated cubic dispersion relation for s{k) follows. It has the property that 
s/ 11^ is a (triple- valued) function of k/ fi only, with parameters a, a^, and /3, 
which set the viscous properties of the fluid and consequently go. 



In the long wavelength limit of the perturbations <^ |//|, the three roots 
behave according to 

s = -2F2^^^ + 0{k), (102) 

s = (/3 + 1)F3A;' + 0{k^), (103) 

s ^ -2iFiijk - ^j^k^ + 0{k^). (104) 

The first mode decays relatively rapidly, as a/^^, on account of the fact that 
F2 is always an increasing function of g^, and hence F2 > 0. The second mode 
is a generalisation of the viscous instabihty (cf. Section 3). For all realistic 
viscous parameters we have -F3(g^) < F2(g^), and with F2(go) = by definition, 
stability of this mode is assured when j3 > —1, the same criterion as for normal 
viscous instability. This mode controls the diffusion of small perturbations 
of wavenumber and amplitude. Lastly, the third mode corresponds to the 
generahsed viscous overstability. It is stable as long as Fi and Fi are both 
positive, which can be easily checked from Eq.(81). In summary, all nonlinear 
wavetrains are stable to small-/c perturbations (in the asymptotic limit of long 
wavelengths that we are examining). 



Numerically, these conclusions can be generahsed to all k. We find that overall 
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stability of the modulational viscous overst ability, and hence the wavetrain, 
can be determined by examining the small- A; limit. In Section 5 in Fig. 6 we 
plot the growth rate of the modulational viscous overstable mode versus k for 
parameters associated with r = 1.5 in Table 1. 



9 Appendix B: Derivation of the Burgers equation 

In this section we provide more details of the derivation of Eq. (45). The proof 
closely follows that of Doelman et al. (2009) who attack a simpler reaction — 
diffusion equation with a methodology developed by Howard and Kopell (1977) 
and Whitham (1974). 



9. 1 Preliminaries 

The nonlinear equation governing the disk, cf. Eqs (14)- (18), we represent in 
a compact vector form: 

dtZi = D,,{a) dlZ, + (Z, d^Z) , (105) 

where Z = (a.u'^.u'y). The diffusion matrix Dij is a function of a only, and 
possesses a first row of zeros (because the continuity equation has no sec- 
ond derivative). The nonlinear vector function depends on both Z and its 
derivative. It comprises advection, rotation, shear, pressure, and the density- 
dependence of the viscosity. 

By defining 6 as earlier, the equation for the nonlinear wavetrain solutions 
Z^{9,^) is 

codeZ^ + fi^ D.^dlZ^^ + Gi (Zo, ^Zo) = 0, (106) 
while the linearised problem for a small perturbation Zi is 

dtZi = udeZi + i^^Dijd^.Zj + iJ.^Dl^d^,Z^ a+^Z,+^i^ deZj. (107) 

Throughout the section we shall refer to the second argument of Gi by the 
vector Also D[ - denotes the a derivative of D^j. Equation (107) is rewritten 
for a single mode oc e^^ so that 

sZi = Cij{e) Zj. (108) 

Having introduced these notational conventions, we next make a few prelimi- 
nary remarks. 
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Note that the null space of the operator dj must include the function deZf. 
That is to say CijdgZj = 0. This is a consequence of the translational symme- 
try of the problem. If Z^{9) is a solution to (106) then so must be Zf{9 + d9). 
Upon substituting the latter into (106), then expanding and linearising in the 
small displacement d9, we arrive at the result. In addition, we assume that 
the null space of Cij is one-dimensional and spanned by dgZf. 

The linear operator possesses periodic coefficients and thus we can make the 
Floquet ansatz: Zj — Zj exp{ik9/ fi). This introduces a new wavenumber pa- 
rameter k upon which the spectrum of Cij smoothly depends. We denote the 
(discrete) eigenvalues of the operator by Sn{k), where n is some index and the 
eigenvalues' functional dependence on k is explicit. Of the various eigenvalues 
associated with the system, we are interested in this section with just one: the 
diffusive mode which governs the relaxation of wavenumber, and which is the 
generalisation of the viscous instability mode. This eigenvalue is zero when 
k = 0, i.e. s(0) = 0, and as we have seen is proportional to k^. It follows that 
its eigenfunction must lie in the null space of Cij when k = 0. Thus we set 
Zj = dgZ^ when k = 0. All the other modes of C^j are assumed to be stable. 

Lastly we denote by Zf^ the nontrivial function in the null space of the adjoint 
of Cij. This we normalise so that 

where we have defined the inner product above by 

(Zi, Z2)= r Z^-Z2d9. 
Jo 

9.2 The modulational equation 

We introduce the slow time scale T and long length scale X of Section 6.2 and 
make the ansatz 

Zi^ z^{9 + e,ii + sex) + s^Zi{9, x, t) 

where the slow modulation Q{X, T) is introduced, and 5 is a small ordering 
parameter. This is substituted into (105) and expanded in orders of 5. The 
result is rather messy and we skip most of the steps. At order 5° we obtain 
Eq. (106), which is satisfied by construction. At order 5^ we obtain the 
derivative of equation (106) which is also satisfied. 

At order 5"^ the balance is 

djZj = QrdeZ^, - A(^) ©xx - Bi{9)Ql, (109) 
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for complicated vectors Ai and Bi. We can be assured this equation is solvable 
if the inner product of the right side with Z'^'^ is identically zero. This furnishes 
the evolution equation for 0, 

Qt - {Z^, A) el - {Z^, B) Qxx = 0. (110) 

What we need to do now is compute the inner products which appear in it. 



The vector may be written as 

+ Di^i^^dee.Z^ + 2f,dlZ^)dy + l^^'D'^d^.Z^ {d^f 
^ 2dZ dZra ^'^^^ ^^^"^ ^ dZ dl^ ^"^^ (M^^m + 9eZ'^) 
+ 2 -^M^ ^^^'^^^ + deZ^)i^deA + deZ'J + ^d^oZ^ 

This expression can be remarkably simplified. By differentiating equation (106) 
twice with respect to /i we obtain an identity which permits us to cancel all 
the terms above save one. We find 

A^-^u;"{fx)doZl 

Therefore, 

{Z^, A) = -ia;"(/.). (Ill) 



The vector may be written as 

dC 

B, = cg d^Z^ + D.jdeZ^ + 2f,D,^d^eZ] + — ^ d.Z^ 

The inner product of this expression with Z^ can also be substantially simpli- 
fied. We turn to equation (106) again and take its first derivative with respect 
to /i; we next make the ansatz Zi = Zjexp{ik6/fi + st), substitute this into 
(107), and take its first and second derivative with respect to k while setting 
k — 0. These three equations yield three solvability conditions from which we 
derive 

||(0) = -2(Z-i, B). (112) 
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In the process of the proof we have needed 

{Z^\ d^Z^) = (113) 

and the identification dkZi = d^Z^ at /c = 0. These proceed from the transla- 
tional invariance of the system. 

The final equation for © is 

Ot + \uj"{ii) Q\ + \s"{Q) Qxx = 0, 

from which we can derive the Burgers equation, Eq. (45). Also, the previous 
appendix supplies us with expressions for a;"(/x) and s"(0) to leading order in 
small II, and this allows us to connect the Burgers equation to the material 
parameters of the disk. Prom Eqs (94) and (103): 

u"{ii) = 2Fi(go), s"{Q) = 203 + l)F,{qo), 

where Qq is the root of F2{q) — 0. 
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